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Experiments have shown that the internal vibrational state of a molecule can affect the intensity 
of high harmonic light generated from that molecule. This paper presents a model which explains 
this modulation in terms of interference between different vibrational states occurring during the 
high harmonic process. In addition, a semiclassical model of the continuum electron propagation is 
developed which connects with rigorous treatments of the electron-ion scattering. 

PACS numbers: 



INTRODUCTION 

In the usual three step model of high harmonic gener- 
ation, HHG is treated as a purely electronic process. A 
single active electron tunnels free of a parent molecule 
and eventually scatters from it, but the molecule itself is 
treated in essentially the same way as a lone atom; no 
more than a complicated potential influencing the active 
electron. However, molecules differ from atoms in an es- 
sential way because they possess non electronic internal 
degrees of freedom, which can themselves be affected by 
the high harmonic process. 

The possibility that such internal degrees of freedom 
could play a detectable role in HHG is intriguing, because 
the intrinsic timescale in HHG - the time necessary to 
ionize, propagate and rescatter - is only half a laser cycle. 
This is faster than many chemically interesting processes, 
holding out the possibility that HHG could serve as a 
probe of molecular motion. Alternatively, tailoring the 
state of a molecule prior to HHG could serve to give 
additional control over the generated light. 

These issues were brought to the fore by an experi- 
ment at JILA 1]. In the experiment, a high harmonic 
generating laser pulse was preceded by a weaker pulse 
whose effect was to stimulate Raman-active vibrations in 
SF 6 molecules. Varying the delay between the two pulses 
was observed to modulate the intensity of the HHG light 
generated by the second pulse. Moreover, the modula- 
tion corresponds to the frequencies of the Raman-active 
normal modes stimulated by the first pulse. None of the 
3 non-Raman-active vibrational modes of SF@ were de- 
tected in the modulated signal. 

This paper presents a quantum mechanical model of 
high harmonic generation in molecules. This model pro- 
vides a framework to interpret the observed modulation 
of high harmonic intensities observed in the JILA exper- 
iment, and is easily extended to systems with more com- 
plicated dynamics. Secondly, it presents a version of the 
three step model which has been improved for the pur- 
pose of treating HHG in molecules with relevant internal 
degrees of freedom. Finally, the modulations predicted 
by this improved model are compared with the modula- 



tions observed in the JILA experiment. This paper reca- 
pitulates and extends work which originally appeared in 



THE VIBRATIONAL WAVEFUNCTION OF THE 
MOLECULE 

An important difference between atomic and molecu- 
lar systems, probed by the aforementioned experiment [l| 
is the presence of vibrational degrees of freedom in the 
latter. For an M-atom molecule, this corresponds to 
N = 3M — 6 (N = 3M — 5 for linear molecules) inter- 
nal degrees of freedom, which can be expressed in normal 
mode coordinates. The vibrational wavefunction of the 
molecule can then be expanded as the product of simple 
harmonic oscillator basis functions in each of the normal 
modes of the molecule: 
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where |niri2 ■ ■ -tin) is the outer product of simple har- 
monic oscillator state \n\) in the first normal mode, \n%) 
in the second normal mode, and so on. 

For the purposes of this paper, all operators will be 
expanded to first order in raising a^W and lowering a' 1 ) 
operators in each normal mode i. At this level of approxi- 
mation, the evolution of the vibrational wavefunction be- 
comes separable, and the overall coefficient A nin2 _. nN (t) 
can be factored into the product of individual coefficients 
a ni (t): 

(t) = a ni (t)a n2 (t) . . .a nN (t). (3) 

(Note that the coefficients a ni (t) should not be confused 
with lowering operators a^.) 

In the interests of simplicity and clarity, the remain- 
der of this paper will use a 1-D picture, describing the 
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evolution of the vibrational wavefunction for each single 
normal mode. The concepts from the 1-D model extend 
readily to the higher intensity regime of coupled modes, 
in the event that operators involving two or more raising 
or lowering operators become important. (23j 

VIBRATIONAL INTERFERENCE 

Although high harmonic generation is primarily an 
electronic process, the vibrational state of the molecule 
can affect the harmonic intensity. This occurs because 
at several points in the high harmonic process, the vi- 
brational wavefunction has amplitudes either to stay un- 
changed or to "hop" up or down from simple harmonic 
oscillator state \n) to states |n± 1), much in the same 
way that a photon in a beamsplitter has amplitudes to 
take two or more paths. As in a beamsplitter experi- 
ment, two or more indistinguishable pathways interfere 
with one another and modulate the output signal de- 
tectably. The multiple pathways at play in the high har- 
monic process are diagrammed in Figure [TJ 

Raman Excitation The first opportunity to change 
vibrational states occurs when the vibrationally cold 
molecules are subjected to the weak, off-resonant initial 
pulse. This causes the molecule to undergo stimulated 
Raman scattering. At the start of the first pulse the 
molecule is in the |0) vibrational state, and it is then 
driven into a coherent superposition of the zeroth and 
first vibrational states 

h^ib) = a (0)|0)+ai(0)|l) (4) 

by the end of the pulse. (For the pulse length and in- 
tensity used in the JILA experiment, calculations show 
no appreciable population of the \2) or higher states af- 
ter the Raman pulse. Accordingly, the |2) and higher 
states have been dropped from this analysis.) Because 
only one normal mode is used in this analysis, the N 
quantum number has been dropped, so that ao = ai o (0) 
and a\ — ai 1 (0). t = is chosen at some time after 
the end of the first pulse, when the stimulated Raman 
scattering is over. 

The vibrational state coefficients follow equations of 
motion given by 



ia ni {t) = 




[QAB«n, + diOtABWrii + la„ i + i + V^i a ",-l)] ■ 

Here ioi is the normal mode frequency, indices A and B 
run over {x, y, z}, E^{t) is the component of the electric 
field in the (body-frame) A direction at time t, Qi is the 
normalized displacement associated with normal mode 
i and aAB(Qi,Q2, ■■■) is the polarizability tensor of the 
molecule. These equations of motion have off-diagonal 
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FIG. 1: The vibrational interference model @] in one dimen- 
sion. The molecule ends the first (Raman) pulse in a super- 
position of the v = and v = 1 vibrational states. After a 
time delay, the two vibrational states are mixed by stimulated 
Raman scattering (transfer matrix M) , "hopping" during ion- 
ization (I) and recombination (R), as well as evolution of the 
ionic wavefunction while the electron is away (N). Interfer- 
ence between adjacent vibrational states modulates the high 
harmonic signal. 



elements only if diCCAB = (2muJi) daAB/dQi\Q l=0 ^ 
0, which is the condition for a mode to be Raman active. 
The polarizability tensor and its derivatives are found by 
performing an unrestricted Hartree-Fock calculation [3| 
using the aug-cc-pVTZ basis set. Q 

Between the two pulses, the |0) and |1) states evolve 
as eigenstates of the simple harmonic oscillator Hamilto- 
nian. The states are then mixed once again by stimulated 
Raman scattering during the high harmonic generating 
pulse. These effects are approximated by a unitary 2x2 
transfer matrix M, where My is the amplitude to be in 
state \i) at the instant of ionization after beginning the 
second pulse in state \j). 

Ionization and Recombination The vibrational wave- 
function evolves further during each of the three steps - 
ionization, propagation, and recombination-of the three 
step model. The molecule's vibrational state hops up or 
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down a level during ionization, evolves while the propa- 
gating electron is away from the molecule, and hops once 
again when the electron recombines with the parent ion. 

These hopping amplitudes arise because ionization and 
recombination, commonly thought of as purely electronic 
processes, are both strongly modulated by molecular dis- 
tortions. This is the simplest way in which the internal 
degrees of freedom in molecules allows for behavior that 
has no analogue in atomic systems. Nonzero derivatives 
of ionization and recombination amplitudes translate di- 
rectly into amplitudes for the molecule to change its vi- 
brational state during these processes. Taylor-expanding 
the ionization operator about the equilibrium configura- 
tion of the neutral molecule, 

df 

/ = / |cq+7^Q + °(Q 2 )> ( 6 ) 



using the identity Q = (a + at )/ V2mw and substituting 

Iq = i\cq, h — (2?7iw) _1 ^ 2 ^, the ionization operator 
can be rewritten 

7 = / + /i(a + at). (7) 

Identical logic gives the recombination dipole vector op- 
erator 

R= R + fli(a + a t ) (8) 

Vibrational Dynamics of the Parent Ion Between the 
times of ionization and recombination, the evolution of 
the internal state of SFjj" is quite complicated. This is 
because SF6 has three degenerate orbitals at the point of 
maximum symmetry. Thus, at any nuclear configuration 
near this maximum symmetry point, these three SF^~ 
orbitals are very nearly degenerate, and are mixed with 
one another strongly by molecular distortions. Orbital 
degeneracies can be broken and orbital energies can cross 
even with relatively small distortions. Because of this, it 
is necessary to treat this interplay between electronic and 
vibrational states when describing the dynamics of SFg~ 
in the vicinity of the maximum symmetry point. 

In its maximum symmetry configuration, SF + belongs 
to the Oh point group, with three degenerate Ti g orbitals 
which transform like axial vectors x, y and z. When the 
molecule is distorted away from the maximum symmetry 
point via either an E g or a T2 g distortion, the triple de- 
generacy breaks up into three nondegenerate electronic 
orbitals. The fully symmetric Ai g or "breathing" mode 
preserves the triple degeneracy. 

The full vibronic (vibrational-electronic) Jahn- Teller 
coupling matrix for a triply degenerate system is given 

byiaa 

(91 - ge + %/% gc, g v \ 

h t = I ffc 9i- 90 - VSge 9i I (9) 

V g v gt 9i + 2 geJ 



This matrix represents the coupling between the states 
with x, y, and z symmetry, caused by vibrational op- 
erators, so that represents the off-diagonal coupling 
between the states with x and y symmetry, g v the off- 
diagonal coupling between the states with x and z sym- 
metry, and represents the off-diagonal coupling be- 
tween the states with y and z symmetry. 

Defining the E g normal mode coordinates Qg and Q e , 
which transform like 2z 2 — x 2 — y 2 and x 2 — y 2 respectively, 
and the T2 g normal mode coordinates Qr) and Qq 
respectively as the coordinates that transform like yz,xz, 
and xy, the functions <?i are given by H 

9i= \K E {Ql + Q 2 e ) + \K T (Ql + Ql + Ql) (10) 
ge = ^V Eg Q e + N E (Q 2 - Q 2 ) + N 1 (2Q 2 - Qf - <$)L) 
ge = ^VE g Qe + 2N E Q e Q e + V^N^Q 2 - Q 2 ) (12) 
Vt 2s Q 6 + NtQ v Qc + N 2 Qz(s/3Q e - Q e ) (13) 

9v = v T 2g Qv + NtQcQz + N 2 Q v {~V3Q t - Qe) (14) 
9C= V T2g Q C + N T QsQ n + 2N 2 Q c Qg (15) 

These constants were found by performing a CASSCF 
state-averaged calculation for the three lowest energy 
states of SF^" for various displacements of the molecule 
away from the maximum symmetry configuration. These 
energies are shown in Figure [H The CASSCF calcu- 
lations were carried using a basis of Hartree-Fock or- 
bitals calculated for neutral SFq. These three adia- 
batic energies were then fitted to the eigenvalues of 
the diabatic Jahn- Teller coupling matrix. This pro- 
cess yielded F T2g = -001209 H/bohr, V" Bg =.1406 H/bohr, 
iVi=-.0362 H/bohr 2 , K T2g = .7288 H/bohr 2 , ^=1.8486 
H/bohr 2 . For the Ai g mode, which does not enter 
into the vibronic Hamiltonian, an adiabatic potential 
E = V Alg Q Alg + \K Alg Q 2 Aig , with V^ le =.0645 H/bohr, 
Ka-i =2.98 H/bohr 2 gives the potential energy surface 
for all three electronic states. 

Here there is a significant distinction between distor- 
tions of type Eg , which break the triple degeneracy but 
have no linear off-diagonal terms, and distortions of type 
T2 9 , which do contribute to off-diagonal coupling. For 
small distortions, an adiabatic electronic state of an E g 
distorted molecule will have the same symmetry - x, y, 
z - as the diabatic electronic states. The adiabatic elec- 
tronic states of a T2 g -distorted molecule, on the other 
hand, are linear combinations of the diabatic orbitals. 
An important simplification is that Vr 2g , controlling off- 
diagonal coupling between different electronic states, is 
small in SF^" and can be neglected for the short time 
between ionization and recombination. 

The evolution of the ionic wavefunction is calculated 
in the vibrational basis of the neutral molecule, but using 
potential energy surfaces calculated for the ion. Poten- 
tial energy curves are found to quadratic order in Q using 
quantum chemistry calculations, then expressed in terms 



4 




A . displacement (bohr) 





\ 


1 1 1 


1 1 / 


993.338 




E=E (] +O.L406J q + 0.S3ISq 2 

E=E (f 0.14061 q + O^ISq 1 

E=E (f 7.3653c- 13 q+ 0.1 1707 q 3 




993.340 








993.342 








993.344 








993.346 




1 


1 , \ 



-0.04 -0.02 0.02 0.04 

E distortion (bohr) 



of raising and lowering operators by substituting Q = 
(2rnw)- 1 /2( a + at), Q 2 = (2mu)- 1 {a + at) (a + at). All 
terms up to linear in raising and lowering operators are 
then used to integrate the time-dependent Schrodinger 
equation to find a transfer matrix N_ describing the evo- 
lution of the ionic wavefunction between ionization and 
recombination. 

Modulation of Harmonic Intensity In the two-state 
model used here, the i— th vibrational wavefunction of 
the neutral molecule after recombination has occurred is 
IVVib) = d |0) + dx where 

(d di)=(oo(0) ai {0)e-^ T )M T L T N T R T . (16) 

Here, e.g. M T denotes the transpose of matrix M . 

The number of photons emitted in a given harmonic 
is proportional to do ■ dg + d\ ■ d\. The high harmonic 
intensity is a sum over all Raman active modes i: 

P{r) = P Q + SiP x (i) cos (lu.t + St) (17) 

The static Pq primarily results from terms of the 
form ao(0)*eto(0), while Pi results from terms of the 
form a (0)ai(0)V WT and ai(0)a o (0)*e- MT . Defining 
W = Mt/tTvt^t _ j} NIM t p Q = ao (o)*Wooa Q (0) and 
Pi cos (ujt + 8) = \{ale lulT Wioao(ti) + c.c). Since h and 
Ri are small relative to I and R , only their first-order 
terms are kept. 

DESCRIBING THE CONTINUUM ELECTRON 
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FIG. 2: a)Spherically symmetric A\ g (breathing mode) distor- 
tions change electronic state energies, but preserve the triple 
degeneracy. Non symmetric E 9 (b) and T20 (c) distortions 
break the triple degeneracy of SFg" at the maximum sym- 
metry point. Adiabatic energies are fitted to the eigenvalues 
of the vibronic coupling matrix (equation |5J) to solve for the 
vibronic coupling constants. 



The evolution of the continuum electron wavefunction 
is strongly influenced by interactions with both the par- 
ent ion and the driving laser. When the electron first 
tunnels free out of the parent molecule, its wavefunction 
is determined by both the molecular potential and the 
electric field of the laser. Once free, it propagates in the 
time-varying field of the laser while feeling a weak force 
due to Coulomb attraction to the parent ion. Finally it 
recollides with the parent ion, and is once again strongly 
distorted by the molecular potential. 

A full solution of the time dependent Schrodinger equa- 
tion for this process would be computationally demand- 
ing for complicated molecules such as SFg. In addition, 
much of the information in the continuum wavefunction 
is not relevant to the HHG problem: only a small part 
of the wavefunction overlaps with the unoccupied orbital 
into which the rescattering electron recombincs. 

One frequently used treatment which avoids the com- 
plications of the full time-dependent Schrodinger equa- 
tion 0, Q is based on a classical or semiclassical prop- 
agation of the continuum electron, ignoring the ionic 
Coulomb potential. The returning electron wavefunction 
is the approximated as a plane wave throughout the re- 
combination process. This approach has been successful 
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in describing the high harmonic cutoff, the chirp of the 
emitted high harmoni c lig ht and other quantities of inter- 



est in atomic systems [lOL UJj. However, as was shown in 
, the plane wave approximation is not adequate to de- 
scribe the returning electron because of the tremendous 
distortion caused by the electron's interaction with the 
ionic potential and by exchange effects with the other 
electrons in the molecule, 13|, [l4( For the time-reversed 
problem of photoionization, it is known that the plane 
wave approximation is prone to error for photoelectron 
energies smaller than the deepest K-shell binding energy. 
Energies attained in high harmonic generation experi- 
ments usually fall below this range. 

This section gives a semiclassical model of the 
free electron propagation which improves upon the 
Corkum/Lewenstcin model by connecting with calcu- 
lated short-range wavefunctions. Any method could be 
used to calculate these short-range wavefunctions. In this 
treatment, the tunneling wavefunction is modeled semi- 
classically using ideas based on the initial value represen- 
tation 15, 3]. The continuum wavepacket of the recol- 
liding electron is described in terms of electron-molecule 
scattering states as used in 12[ , calculated in the absence 
of an external electric field. 

In this way, the molecular potential nontrivially affects 
the electron wavefunction at the two times when the elec- 
tron is near the molecular ion. When the electron is far 
from the molecule, the comparatively simple evolution 
of its wavefunction is described using the shortest-time, 
dominant contribution to the Gutzwiller propagator [fjj • 
Finally, stationary phase arguments serve to identify the 
isolated trajectories that encapsulate the effect of the 
electron propagation in the field on high harmonic gen- 
eration, greatly reducing the computational burden of 
propagating the continuum electron wavefunction. 



Tunneling Ionization 

During the ionization step, the tunneling electron 
wavefunction is described in a simple 1-D WKB tunnel- 
ing picture, in which electrons arc allowed to tunnel only 
in directions parallel to the laser's applied electric field. 
This approach is motivated by the semiclassical "initial 
value representation" [IB], [l6| , where a source wavefunc- 
tion acquires an imaginary phase (and hence an expo- 
nential growth or decay) along a trajectory that passes 
through a classically forbidden region. The unperturbed 
highest occupied molecular orbital (HOMO) here serves 
as the source wavefunction, so that the tunneling wave- 
function is approximated by the unperturbed HOMO in 
the classically-allowed region near the molecule, connect- 
ing to a WKB exponential which is set equal to the 
HOMO at the inner turning point and decays exponen- 
tially until it reaches the outer turning point. SFq has 
three degenerate HOMOs: one of these is illustrated in 



FIG. 3: One of three degenerate orbitals of SF6. Red denotes 
positive lobes; blue denotes negative lobes. 




FIG. 4: The tunneling wavefunction is approximated as an 
unperturbed molecular HOMO inside the classically allowed 
region, connecting to a decreasing WKB exponential in the 
classically forbidden region. Stationary phase trajectories 
leave from the outer turning point, beginning with zero ve- 
locity at time of ionization. 



Figure [3] 

The classically forbidden region, illustrated in Figure 
21 is defined by two turning points, a and 6, between 
which V(x) ~ E > 0. The slope of V(x) is C\ at inner 
turning point a and Ci at outer turning point b. The 
tunneling wavefunction ipt{r, t) is now found by applying 
WKB connection formulas. The appendix derives the 
ratio of the tunneling wavefunction at the outer turning 
point to the tunneling wavefunction at the inner turning 
point 



ip(x = b) 
ip(x = a) 



± P -r|Hi|i/6. Bi (°) 



I — I 



Ai(0) 



(18) 



In this approximation, the tunneling wavefunction be- 
haves like an Airy Bi function near the outer turning 
point, i.e. it has no linear complex phase term. This 
property will be revisited in the section dealing with sta- 
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tionary phase analysis. 



Semiclassical Propagation 

When the active electron has tunneled free from the 
molecule, the evolution of its wavefunction is controlled 
by the oscillating electric field of the laser, plus a residual 
Coulomb attraction to the molecular ion. This relatively 
simple evolution continues until the electron returns to 
the molecule, when the complicated molecular potential 
again becomes significant. During this excursion, un- 
til it re-enters the non-Coulomb part of the potential, 
the continuum wavefunction can be approximated using 
Gutzwiller's semiclassical propagator [l7| : 



K(r,t;r ,to) = ^ (27ri) _3/2 y/C(r, t; f , t ) x 

cl. traj. (19) 

exp(iS{f,t;f 0l t ) - i<f>) 

Here S(r, t; To, to) is the action integral S = 
J L(q,q,t)dt calculated for a classical trajectory start- 
ing at (fo,to) and ending at (r, t), while C(f,t;fo,to) = 
| j drel, where ta is the A-component of the vector 
r. is a phase factor equal to \ times the number of 
conjugate points crossed by the trajectory [17j. 

This propagator acts on the tunneling wavepacket ipt 
to give a semiclassical continuum wavepacket 

i>cif,t) = J d 3 r J dt K{r,t;r ,t )Mro,to)- (20) 

until the molecular potential asserts itself during the 
rescattering. 

During the terminal portion of the scattering process, 
when the scattering wavefunction acquires its maximum 
dipole matrix element with the molecular HOMO into 
which it rccombines, the electronic wavefunction is ex- 
panded into a truncated (but in principle complete) basis 
of field-free electron-molecule scattering orbitals, calcu- 
lated using techniques described in references [l^.ll9l[20|. 
Beyond the range of the molecular potential, i.e. for 
r > ro, the (l,m)-th independent scattering state is ex- 
pressed as a partial wave expansion in terms of incoming 
and outgoing Coulomb radial functions f Ei (r) and the 
scattering S-matrix as 



ipE,lm(r) = -^=f El (r)Yi m (6, 



I'ttv 



(E),r > r . 



(21) 



The laser electric field is typically far smaller when the 
electron returns to the ion than it was when it departed. 
Even if this were not the case, this external force is less 
than that due to the electron-ion interaction force when 



the electron is in the molecular field. Neglecting the ef- 
fect of the external field on the electron during its brief 
recollision with the ion, the time-dependent wavefunction 
becomes 

Vw(r,t) = J dEY,Ai m (E)^EMn(r)e- lEt (22) 

lm 

and the expansion coefficients Ai jTn (E) are given by 

A lm (E) = e lEt ( d 3 r<p*E, lm (r>c(r, *) (23) 



for some chosen time t when ip c is projected onto the 
scattering states. This projection time is chosen so that 
the bulk of the returning wavepacket is approaching close 
to the ion but has not yet recollided. 
The dipole recombination amplitudes 



dl, m (E) = (lpB,lm\ * ' ^IV'HOMO) 



(24) 



are calculated between the distorted scattering states and 
the molecular HOMO into which the electron recombincs. 

Equations [5D] and [55] now define the expansion coeffi- 
cients Ai m (E) for some chosen projection time t, in terms 
of two three-dimensional integrals over initial and final 
positions, one integral over the time of ionization, and 
a summation over all possible classical trajectories. Sta- 
tionary phase arguments dramatically simplify the calcu- 
lation of these expansion coefficients. 



Stationary phase calculation of scattering coefficients 

A difficulty which must be resolved when using the 
Gutzwiller propagator to describe the evolution of the 
continuum wavefunction far from the molecule and scat- 
tering/tunneling wavefunctions to describe its evolution 
near the molecule is the very different physical pictures 
employed in the two treatments. The Gutzwiller prop- 
agator K(r,t;fo,to) defined in equation [111] involves a 
summation over the classical paths having (Fo,£o) anc ^ 
(r, t) as their endpoints, whereas short-range treatments 
of scattering and tunneling use more familiar wavefunc- 
tion descriptions. Moreover, previous semiclassical treat- 
ments such as [8, 9] have found only a few such classical 
paths to be significant in describing the high harmonic 
process. 

Both of these difficulties may be resolved by using sta- 
tionary phase techniques to search for families of classical 
trajectories which add together to give nonvanishing con- 
tributions to the expansion coefficients Ai_ rn {E) defined 
in equation I23I 

Consider the contribution to Ai^ m (E) made by all clas- 
sical paths originating at a given point (ro,to)- Making 
the usual assumption of isolated trajectories, there will 
be one classical path connecting starting point (ro>^o) to 
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any given ending point (r , t) , which will have the prop- 
erty of extremizing the action integral S(f,t;fo,to) = 
J L(q, q, t)dt with respect to any perturbation which does 
not alter the starting or finishing points. The connection 
to the wavefunction treatment of the scattering and tun- 
neling wavefunctions is made by noting that S(f, t; r , to) 
appears as a phase term in the Gutzwiller propagator. 
Thus, the change in accumulated action which is caused 
by changing the endpoints of a classical path corresponds 
to a change in the phase of the propagated wavefunction. 

This variation of the action integral with respect to the 
starting and ending points of a trajectory, is known (see 
Goldstein et al [2l[ section 8.6) as the A variation of the 
action, given by 



AS(r, t;r a ,t ) = {p A Sq/ 



(25) 



The condition for a nonoscillating integrand in Eq. [23] 
is now that the variation of the phase of the Gutzwiller 
propagator resulting from Eq. [25] must be offset by the 
variation of phase of the tunneling or scattering wave- 
function. In this paper, a trajectory where the A varia- 
tion of the action is counterbalanced by the phase of the 
tunneling wavefunction at (fo, to) and by the phase of the 
scattering wavefunction at (r, t) will be known as a sta- 
tionary phase trajectory." The initial and final points 
(fo,to) and (r,i) identify the points in the 7-dimensional 
integral of Eq. [23] where the integrand oscillates slowly, 
giving a non-canceling contribution to the expansion co- 
efficients Ai m (E). 

Stationary phase trajectories will be found to corre- 
spond with the trajectories used in prior semiclassical 
theories of high harmonic generation. However, the cur- 
rent approach allows for more detailed treatments of 
the short-range tunneling and scattering wavefunctions, 
where semiclassical methods may give unsatisfactory de- 
scriptions of the physics. 

Stationary phase trajectories may be found by expand- 
ing the phase-oscillating parts of the 7D integral from Eq. 



A lm (E) = / d 3 r / d 3 r / dt e lEt y/C(r, t; f , t ) 



exp(iS(r,t;r ,t ) - i<P)ipE,i m (r,t)')pt{ i r'o,to) 



(26) 



about the starting point (fo c ,ioc) an< ^ arj out the ending 
point (r c ,t) 

The A variation gives the expansion of the action in- 
tegral 

S(f c + 5r,t; f 0c + Sr 0l t c + St ) = S(f c ,t; r 0c ,io c )+ 
VaOt a + ttt; — — or A or B + PoAdr A+ 



2 drAdrs 
■oroAOroB + Hdt Q + --^-(dto) ■ 



2 dr Adr B 



The condition for a nonoscillatory integrand is now 
that the first-order terms in 5ro,Sf and 6to must disap- 
pear. 

For the tunneling wavefunction, which resembles a de- 
clining WKB exponential in the forbidden region and has 
no oscillatory component, this corresponds to a trajec- 
tory which leaves the molecule with zero initial momen- 
tum. 

For the scattering wavefunction, ignoring the angular 
derivative of the spherical harmonics, the phase evolution 
of the scattering states is given by the asymptotic form 
of the Coulomb wave functions (r) 

r E ,i m if c + Sr) = { f£(r c )Y? m (9 e , <j> c ) exp(zfc(r c )«5r)- 

-^/^^ m ,(^,^)^ m / ;im ( J E;)ex P (-ifc(r c )fr)) 

frZi ~ z v 2 

(28) 



V ,m' 



where ki(r) — \f(2(E — Vi(r))). The condition for a 
nonoscillatory phase integrand is that p r = —ki{r) — 
y/(2(E - Vi(r))). 

Neglecting the angular derivatives of the spherical har- 
monics is justified because the action along the recollision 
trajectory is much greater than that in a low-order an- 
gular solution. It translates into a stationary phase con- 
dition that the trajectory must return with zero angular 
momentum. 

Finally, setting the coefficients of Sto to zero yields the 
condition that 



^ L + V(ro,to) = E KO MO- 
2m 



(29) 



(27) 



Thus, a stationary phase trajectory is launched with 
zero momentum from the classical turning point and re- 
turns to the molecule with zero angular momentum, and 
with kinetic energy equal to the energy of the scattering 
state. This is a familiar result from, e.g., @, Q, with the 
distinction that the present work considers the effect of 
the electron-ion Coulomb interaction during the contin- 
uum propagation of the electron. Also, the wavefunction 
is projected onto scattering states shortly before recolli- 
sion, rather than treating the electron in a Volkov approx- 
imation throughout the recollision with the molecular 
ion. Because the electron-molecule scattering states are 
calculated with no external electric field present, there is 
a slight dependence on the time at which the wavefunc- 
tion is projected onto scattering states; here the projec- 
tion is made when Lot for the laser cycle is equal to 3.9, 
i.e. when the short trajectories with energies equal to the 
energy of the 39th harmonic have nearly returned to the 
molecule. 

Because the linear phase variation of the integrand 
vanishes in the vicinity of a stationary phase trajec- 
tory, the expansion coefficients Ai m {E) are now found 
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via Gaussian integrals: 



A lm (E) = (2m)- 3 / 2 ^=f E ;(r c )Y: m (9 c ^ c )x 



exp(iS(f c ,t;fo c ,t 0c ) - i<j))4> t (roc)e' 



-iEHOMOtoc 



((/) = I d(^ )exp(-^ T (<5t ) 2 ))x (30) 



((II) = /d 3 (^ )exp(- 



2 droAdroB 



Sr 0a 5r 0B ))x 



((///)= [ d 3 (6r)eM % 7;^-6r A 6r B )) 
J 2 dr A dr B 

where the integrals labeled (I), (I I) and are eval- 

uated as 



(II) = (2mf/ 2 \ 



(I) = V2m\^\ 

d 2 s 



dt? 



(Ill) = (2«) 3/2 | 



dr 0A dr 0B 
d 2 S 



dr A dr B 



/2 
1/2 
1/2 



(31) 
(32) 
(33) 



yielding expansion coefficients 



A lm {E) = {2^)\\^y, 
,d 2 S 1/2 j n 

V2' 



(34) 



/Br( r c) y im( 6 'c,</'c)exp(iS'(r c ,t;ro c ,to c ) - i0)x 

V' i (roc)e~ 4£HOMO * 0c 

Once these expansion coefficients have been calculated, 
the dipole matrix element between the distorted scatter- 
ing wave and the molecular HOMO is simply 



(35) 



111) 



COMPARISON WITH EXPERIMENT 

The model of vibrational interference connects with 
this treatment of the high harmonic process when IR is 
set to D(E). This is broken up by setting I — ipt(ro,to) 

and R = D(E) / 1. Both of these quantities are calculated 
for a molecule at the equilibrium geometry, and for a 
molecule displaced by 0.1 bohr in the normal mode coor- 
dinate. This involves recalculating the scattering states 
and recombination dipoles for each distorted molecular 
geometry. 

The modulation of the 39th harmonic was chosen for 
purposes of comparison with experiment, since this har- 
monic was considered in detail in [l[ . The 39th harmonic 





FIG. 5: Population of |0) and |1) vibrational states after the 
high harmonic process as a function of angle for the Ti g nor- 
mal mode transforming like xy. 



falls close to the measured cutoff, and can only be pro- 
duced by a half-cycle coming close to the maximum of 
the gaussian envelope of the laser pulse. 

The JILA experiment used a gas jet as a source of SFg, 
giving no preferred molecular orientation. However, both 
ionization and recombination amplitudes are highly de- 
pendent on orientation. Therefore, a rotational average 
was calculated for both the static and oscillatory parts of 
the harmonic intensity. Only those polarizations perpen- 
dicular to the propagating laser beam for a given molec- 
ular orientation were included in these averages. 

A previous workQ compared only the angular aver- 
aged modulation of the entire signal to experiment, by 
calculating the modulation to first order in raising and 
lowering operators. For the present work, the modula- 
tion at a particular molecular orientation was calculated 




5 5.5 1 1.5 E Z.5 3 



FIG. 6: Modulation of HHG signal resulting in final vibra- 
tional state |0) as a function of angle for the Ti g normal mode 
transforming like xy. a) real component of modulation b) 
imaginary component of modulation 

using the full transfer matrices I, R and M. The re- 
sulting vibrational state populations and modulations of 
the different components of the HHG signal show a rich 
angular structure which is lost upon angular averaging. 

Note that all calculations presented here have used the 
separable approximation of Eq[3] for the amplitudes of 
the different vibrational modes. This is expected to be 
accurate to the extent that the ground vibrational state 
dominates, meaning that for the amplitudes that we have 
calculated here, the higher-order nonseparable pathways 
are at least beginning to become important, which dimin- 
ishes the validity of this approximation. Nevertheless, it 
is still expected to have at least qualitative and perhaps 
even semi-quantitative validity for the range of parame- 
ters studied here. 
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FIG. 7: Modulation of HHG signal resulting in final vibra- 
tional state |1) as a function of angle for the Ti g normal mode 
transforming like xy. a) real component of modulation b) 
imaginary component of modulation 

Using the separable approximation, the final popula- 
tions of the |0) and |1) for each normal mode and their 
modulation as a function of time were calculated as a 
function of angle using equation 1171 Figures I5I6I7I show 
final state populations and the real and imaginary com- 
ponents of the modulation as a function of angle for the 
Tig mode transforming like xy. Figures I8I9I10I show fi- 
nal state populations and modulations for the E g mode 
transforming like 2z 2 — x 2 — y 2 , while Figures [111121131 
show final populations and modulation fractions for the 
totally symmetric A\ g mode. In all three modes, the 
population of the |1) vibrational state is modulated much 
more heavily than the population of the |0) state. 

A noteworthy feature of these figures is that in re- 
gions where Raman excitation during the HHG pulse is 



10 





FIG. 8: Population of |0) and |1) vibrational states after the 
high harmonic process as a function of angle for the E g normal 
mode transforming like 2z 2 — x' 2 — y 2 . 



weak, the modulation of the HHG signal due to d^do (ie, 
HHG processes ending with the molecule in the ground 
vibrational state) tends to cancel modulation due to d\d\ 
(leaving the molecule in the first vibrational state). For 
this reason, angle-averaged modulations of the overall 
HHG signal due to the T 2s and E g normal modes give 
almost zero overall modulation, while the modulations 
of the d^do or d*d% give a large fractional modulation. 
(The Ai g mode experiences strong Raman excitation at 
all molecular orientations, and experiences less cancella- 
tion as a result.) 

It is not clear why this cancellation is not observed in 
the experiment, where the T2 9 mode is typically the most 
visible [22| . This could arise due to a variety of causes, 
such as preferential detection of the d\d\ component of 
the signal relative to the d^d^ component, or some mech- 



FIG. 9: Modulation of HHG signal resulting in final vibra- 
tional state |0) as a function of angle for the E g normal mode 
transforming like 2z' 2 — x' 2 — y 2 . a) real component of modu- 
lation b) imaginary component of modulation 



anism changing the phase between the modulation of the 
two components and thereby eliminate the cancellation. 
Table Mil shows the fraction of molecules finishing the 
HHG process in the |0) and |1) states, while table Q] com- 
pares the modulation of the d^do and d\d\ components 
of the signal, the modulation of the total signal, and the 
two experimental runs which were able to detect modula- 
tions at all three vibrational frequencies. Tables [IV] and 
HIl shows the same information when the HHG process 
is allowed to populate vibrational states up to |4) . It is 
apparent from Table HVl that including higher vibrational 
states of the T 2g and E s modes do not greatly affect the 
calculated modulations. For the A\ g mode, which under- 
goes stronger Raman excitation, the presence of \2) and 
higher vibrational states does become important. 
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FIG. 10: Modulation of HHG signal resulting in final vibra- 
tional state |1) as a function of angle for the E g normal mode 
transforming like 2z 2 — x 2 — y 2 . a) real component of modu- 
lation b) imaginary component of modulation. 



FIG. 11: Population of |0) and |1) vibrational states after the 
high harmonic process as a function of angle for the totally 
symmetric Ai g mode. 



TABLE I: Peak-to-peak Modulation, Theory vs. Experiment 
(2 state model) 



Mode Experiment 1 Experiment 2 


|0> 


|1> 


Total Signal 


Ai s .06 .105 


.0439 


.187 


.0246 


T 2g .105 .122 


.0010 


.236 


.0016 


E s .025 .029 


.0255 


.121 


.0020 



TABLE II: Peak-to-peak Modulation, Theory vs. Experiment 
(5 state model) 



Mode 


10} 


|1> 


12} 


|3> 


14) 


Total modulation 


Ai 9 


.0495 


.249 


.508 


1.01 


1.44 


.0901 


T 2 g 


.0099 


.232 


.708 


1.45 


1.84 


.0011 


E 9 


.0264 


.142 


.907 


1.55 


1.88 


.0071 



Although the agreement with experiment is not per- 
fect, it is nevertheless significant that the simple model 
of vibrational interference presented here agrees with ex- 
periment to the correct order of magnitude. This is par- 
ticularly notable in light of conventional Raman spec- 
troscopy, in which the Ai g peak is 20 times more promi- 



nent than the others. It is difficult to precisely gauge the 
agreement of theory and experiment, due to the paucity 
of experimental data. Peak-to-peak modulations vary ex- 
tensively from one experimental run to another 22| with 
this SFg experiment. The modulation at 525 cm -1 , cor- 
responding to the T2g mode, appears most prominently 
in the experimental data, yet it gives the smallest modu- 



12 




FIG. 12: Modulation of HHG signal resulting in final vibra- 
tional state |0) as a function of angle for the totally symmet- 
ric Ai g mode, a) real component of modulation b) imaginary 
component of modulation 



TABLE III: Vibrational State Population After HHG process 
(2 state model) 

Mode |0) pop. |1) pop. 





0.71 


0.29 


T 29 


0.87 


0.13 


E 9 


0.69 


0.33 



lation in this treatment. The prominence of the T2 9 mode 
modulation may suggest that the off-diagonal Jahn- Teller 
coupling Vr 2g is larger than obtained in the present cal- 
culations. Alternatively, it may be necessary to model 
the experiment in more detail - i.e. to incorporate the 
spatially-varying laser intensity, the uncontrolled carrier 
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FIG. 13: Modulation of HHG signal resulting in final vibra- 
tional state 1 1) as a function of angle for the totally symmet- 
ric Ai g mode, a) real component of modulation b) imaginary 
component of modulation. 



TABLE IV: Vibrational State Population After HHG process 
(5 state model) 



Mode 


0) pop. 


|1) pop. 


|2> pop. 


|3) pop. 


[4) pop. 


Ai s 


.611 


.310 


.0719 


7.03xl0~ 3 


4.37xl0~ 4 


T 29 


.955 


.0438 


6.92xl0" 4 


7.90xl0 -6 


8.73xl0~ 8 


E s 


.826 


.168 


5.91xl0" 3 


1.63xl0" 4 


3.59xl0~ 6 



envelope phase, the combination of multiple laser half cy- 
cles, etc - beyond that which has been included in the 
present theoretical description. 
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CONCLUSIONS 

The problem of high harmonic generation in molecules 
can be conceptually separated into two parts: the evo- 
lution of the continuum electron, and the evolution of 
the internal (vibrational) wavefunction of the parent ion. 
This paper describes the evolution of the continuum elec- 
tron in a model which combines a semiclassical treatment 
of the propagation with a fully quantum mechanical de- 
scription of the electron-molecule scattering. This flex- 
ible and robust model has a simple conceptual link to 
existing semiclassical models, yet it allows for a sophis- 
ticated treatment of the complicated electron-molecule 
scattering. The internal dynamics of the parent ion are 
tracked throughout the high harmonic process. Together, 
these two innovations serve to give an unprecedented 
view of high harmonic generation in a comparatively 
large, complicated molecule with many internal degrees 
of freedom, giving results which agree with experiment 
to within an order of magnitude. 

The possibility that high harmonic generation may 
serve as an ultrafast interferometric probe of a molecular 
vibrational wavefunction is extremely promising. Such a 
wavefunction need not be prepared by an initial Raman 
pulse, as was the case for the JILA experiment. Instead, 
a preparatory pulse could photoionize a molecule, excite 
it to a higher electronic state, or trigger the beginning 
of some other chemical process. In this way, vibrational 
wavepacket evolution during chemical processes could be 
observed as it happens. 
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APPENDIX: TUNNELING IONIZATION 

In this approximate treatment, the wavefunction in 
the forbidden region is found using the WKB connec- 
tion formulas. As is illustrated in Figure 21 the classi- 
cally forbidden region is defined by two turning points, a 
and b, aligned in the downfield direction, between which 
V(x) -E>0. The slope of V{x) is C x > at the inner 
turning point a and Ci < at the outer turning point b. 

Near turning point a, k 2 (x) = 2m{E — V{x)) w C\{a — 
x) and the time independent Schrodinger equation is 

ijj" + k 2 {x)4> = 0, (36) 

which has solutions near x = a of 

V>(s)=Ai( \ 2/3 > i+Bi( > 2 (37) 



where Ai and Bi are Airy functions, asymptotically be- 
having like 



(2 7 r)- 1 /2 z -i/4 exp [^ z 3/2 ] (38) 

Ai(z) =>^_oo ^- 1 /2(„ z )-l/4 sin( 2 ( _ z)3 /2 +7r/4l39) 
7T- 1 /2^-V4 exp ( 2^3/2) (4Q) 
_«, ^- 1 /2(- Z )-l/4 cos (| ( _ 2 )3/2 + ^4!) 



M(z) 
Vi(z) = 

Bi(z) 
Bi(z) = 



where z = C\^{x — a) 
Similarly, near x = b 

ip(x) = Ai(C 2 1/3 (6 - x))di + Bi(C 2 1/3 (6 - x))d 2 (42) 

Under the barrier but away from the turning points, 
the WKB wavefunction is given by 

$(x) = Tr~ 1/2 \k{x)\- 1/2 exp( / \k{x')\dx') sin(0)- 



l\k(x)\-^ 2 exp(- / \k{x')\dx') cos(0) 



(43) 

for some value of <f>. 

Setting \k{x'\dx' = Y, note that 

\k(x')\dx' = T - [ \k{x'\dx' (44) 

J X 

and for x close to b 

k(x')dx' = -Cl /2 {x-bf' 2 



(45) 



Connecting the asymptotic forms of the Airy functions 
to the WKB solution in the forbidden region gives the 
solution for x ~ b as 

^(x) = xxb 2C 2 ~ 1/6 sin(0)e r Ai(C 2 1/3 (fe-a;))- 
ic 2 - 1/6 cos(0) e - r Bi(C 2 1/3 (6-.x)), 

giving a scattering phaseshift <5 of 

S = tan~ 1 (4e 2T tan(</>)) 

yielding a resonance at <f> m 0, S w 

Similar logic gives the wavefunction for x ~ a 

i>{x) = XRSQ C- 1/6 sin(</,)Bi(C 1 1/3 ( : r - a))- 



(46) 



(47) 



C 1 1/6 cos(0)Ai(C 1 1/3 (a;-a)) 



(48) 



Finally, setting (j) = 0, the ratio of the tunneling wave- 
function at the outer turning point to the tunneling wave- 
function at the inner turning point is 

« = k r <'!" I,e S < 49 » 

ip{x = a) 2 62 Ai(0) 
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